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ABSTRACT 


Theory  for  wave  propagation  in  inhomogeneous  media  is  used  in 
a critical  examination  of  observed  P-wave  travel  time  and 
amplitude  anomalies.  The  NORSAR  array  data  analysed  exhibits 
several  features  predicted  theoretically  assuming  the  crust 
and  upper  mantle  to  be  randomly  inhomogeneous.  It  is  demonstrated 
that  in  modelling  plane  wave  time  deviations  a better  fit  to 
the  observations  is  obtained  by  including  a stochastical  term 
in  the  plane  wave  equation.  It  is  also  demonstrated  that  wave 
scattering  effects  is  of  critical  importance  for  seismic  velocity 
estimation  when  using  data  from  small  aperture  arrays.  In  such 
cases,  a spatial  sampling  of  the  area  under  investigation 
is  recommended  in  order  to  obtain  reliable  velocity  estimates. 
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1 . INTRODUCTION 

• 1 Anomalies  --  Deviations  from  Uniform  Plane  Seismic  Waves 

Signal  similarity  across  seismic  arrays  like  NORSAR  is  an 
important  theoretical  assumption  which  fails  to  be  strictly 
valid  in  practice.  This  report  deals  with  the  deviations 
from  uniformity  encountered  at  this  seismological  observatory, 
and  tries  to  explain  them  by  random  irregularities  (inhomo- 
geneities) in  the  crust  and  upper  mantle.  The  energy  received 
at  NORSAR  seismometers  from  events  at  epicentral  distances 
greater  than  30  is  usually  treated  in  terms  of  uniform 
plane  waves.  For  n equivalent  sensors  comprising  the  array 

we  expect  the  waveform  to  follow  the  mathematical  formulated 
law: 


Ai=Ao  i = 1 , . . . , n (1.1.1) 

Ti  = ri*u  i = 1, . . . ,n  (1.1.2) 

where  A.,  is  the  signal  amplitude  at  the  i-th  sensor,  h is 
the  uniformly  expected  amplitude,  xi  is  the  time-shift  at 
the  i-th  sensor  relative  to  a reference  time,  r.  is  the 
position  vector  of  the  i-th  sensor  relative  to  a reference 
point,  and  u is  the  travel  time  gradient  (slowness) . 

In  this  report  words  like  anomaly,  deviation,  residual, 
fluctuation  and  error  will  be  used  synonymously  to  express 
the  imperfectness  of  the  observed  wavefield  across  NORSAR 
relative  to  a uniform  plane  wave  described  by  eq,  (1.1.1) 
and  eq.  (1.1.2).  Fig.  1.1  shows  a small-scale  picture  of 
the  sensor  configuration  within  subarrays  04B,  05B  and 
08C  at  NORSAR. 


Preceding  page  blank 
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Fig.  1.1  Section  showing  the  configuration  of  single  sensors  within 
NORSAR  subarrays.  Each  black  dot  gives  the  position  of  a 
vertical  short-period  seismometer. 

The  whole  array  consists  of  22  suoh  subarrays  (see  Fig.  4.3) . 
Fig.  1.2  illustrated  what  kind  of  deviations  actually  are 
observed  at  NORSAR.  Delaying  traces  according  to  plane  wave 
does  not  give  a proper  line-up.  The  amplitudes  also  show  a 
non-uniform  picture.  Both  these  facts  are  in  disagreement 
with  the  uniform  plane  wave  model. 

It  should  be  stressed  that  the  time  anomalies  on  subarray 
level  c're  physically  real,  not  purely  measurement  error, 
and  that  amplitude  differences  are  too  big  to  be  attributable 
to  geometrical  spreading  or  elastic  attenuation. 
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Fig.  1.2  Delayed  single  instrument  recordings  for  an  Aleutian  event 
plotted  ir.  the  same  frame  according  to  the  best  plane  wave 
fit,  denoted  PL.W,  and  the  more  correct  line-up  determined 
by  correlation,  denoted  C.  Sampling  rate  is  20  Hz. 


1 . 2 Travel  Time  Anomalies 

Fig.  1.3  shows  the  distribution  of  time  differences  between 
a 1 ine-up  determined  by  the  best-fitting  plane  wave,  and  the 
"correct"  line-up  determined  by  correlation.  The  distribution 
seems  to  have  a Gaussian  shape.  It  should  be  noted  that  this 
distribution  was  obtained  by  considering  the  plane  wave  front 
deviations  within  subarrays  and  thus  representative  only 
for  deviations  from  the  plane  wave  shape  within  a small 
portion  (7-8  km)  of  the  wave  front. 

The  pattern  of  subarray  time  anomalies  shows  no  systematic 
trend  as  function  of  epicenter  region.  This  is  demonstrated 


xjc  -nn  -ix  o kb  /oo 

pl.u.residuhl  tnstC) 


IX 


u 


1*3 


Hx 


Jo 


Fig.  1.3  Distribution  of  plane  wave  residuals  on  subarray  level.  The 
correct  line-up  is  determined  by  the  correlation  procedure 
described  by  Gangi  and  Fairborn  (1968) . Well-recorded  events 
filtered  by  a Butterworth  bandpass  filter  in  the  range  0.8-2. 8 Hz 
have  been  used. 


by  the  random  distribution  of  plane  wave  residuals  of  sub- 
arrays 05B  and  12C  in  Fig.  1.4.  On  the  other  hand,  the 
residuals  are  stationary  in  the  sense  that  they  repeat  them- 
selves for  events  from  the  same  epicentrul  area.  This  holds 
even  for  events  from  regions  of  great  structural  complexity 
(e.q.,  continental  margins  of  the  Pacific  Ocean)  as  demon- 
strated for  subarray  03C  in  Fig.  1.4  The  stationarity  seems 
to  indicate  a near-receiver  structure  causing  the  anomalies. 
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Fig.  1.4  Plane  wave  residuals  of  the  tyoo  whose  distribution  is 

05BWandni2n  ‘ 3'  demonstrating  the  randomness  (subarrays 

® " d 12  } and  the  stationarity  (subarray03C)  in  the 

P« ttprn . 


1 • 3 Amplitude  Anomalies 

The  amplitude  pattern  across  NORSAR  has  been  subject  to 
a statistical  analysis  (Husebye  et  al,  1974),  where  one  of 
the  important  conclusions  turns  out  to  be  the  stationarity 
of  the  amplitude  pattern  for  events  from  the  same  epicenter 
region.  This  observation  agrees  with  the  travel  time  observa- 
tions. As  for  the  time  residuals,  we  are  inclined  to  favor 
near-receiver  structures  causing  this  effect. 

The  distribution  of  maximum  signal  amplitudes  across  the 
132  equivalent  short  period  vertical  seismometers  at  NORSAR 
is  skewed  to  the  right,  i.e.,  higher  values.  Statistical 
tests  favored  a lognormal  distribution,  which  was  explained 
by  Ringdal  et  al  (1972)  as  the  multiplicative  effect  of 
a multilayered  earth  structure.  Tatarski  (1961)  explains 
a lognormal  amplitude  distribution  in  terms  of  wave  propa- 
gation in  a random  mecium. 
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1*4  The  Choice  of  a Random  Medium  Model 

Berteussen  (1975)  discusses  the  possibility  of  explaininq 
time  delay  corrections  at  NORSAR  (i.e.,  anomalies)  by 
interfaces  in  the  crust  and  upper  mantle.  He  concludes  that 
current  deterministic  crust  and  upper  mantle  models  are  able 
to  explain  only  a small  part  of  the  observed  anomalies.  In 
a search  for  a better  approach  to  the  problem  of  explaininq 
anomalies,  we  have  found  it  reasonable  to  consider  a random 
medium  model  as  described  theoretically  by  Chernov  (1960)  . 

This  approach  has  been  widely  used  in  other  fields  of  physics, 
especially  in  acoustics  and  radio  telemetry.  The  choice  of  a 
random  medium  model  is  supported  by  the  existence  of  the 
Gaussian  travel  time  anomaly  distribution  and  the  lognormal 
amplitude  distribution  previously  mentioned. 

This  report  will  contain  a condensed  theoretical  presenta- 
tion of  P-wave  scattering  in  a randomly  inhomogeneous  medium. 
The  applicability  of  the  model  for  the  crust  and  upper 
mantle  beneath  NORSAR  will  be  discussed  and  compared  to  simi- 
lar results  obtained  by  Aki  (1973)  and  Capon  (1974)  for  the 
Large  Aperture  Seismic  Array  (LASA)  in  Montana,  and  also  by 
Capon  and  Berteussen  (1974)  for  NORSAR  data.  Also  a refined 
version  of  the  plane  wave  model  (eq.  (1.1.1)  and  eq.  (1.1.2)) 
proposed,  accounting  for  the  stochastic  part  of  wave 
field  caused  by  random  scattering. 
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2 • THE  RANDOM  MEDIUM  CONCEPT 
2 • 1 Random  Inhomoqene i tie s 

Deviations  from  homogeneity  occur  in  every  real  medium. 

Usually  observed  inhomogenei Mes  are  classified  in  two  main 

types,  regular  or  random.  A spatial  mapping  of  a real  stationary 
medium  could  be  briefly  formulated  as 


a<r>  = a0  + al(?)  + 


(2.1.1) 


where  a is  some  physical  parameter  describing  the  medium, 
aQ  is  a constant,  is  the  regularly  varying  term,  and  E 
is  the  unknown  random  deviation  from  the  deterministic  part 

a0  + a1(r).  If  E denotes  expectation,  e(r)  is  required  to 
satisfy  the  condition 


E{e(?) } = 0 


(2.1.2) 


If  regular  changes  in  a real  medium  occur  in  some  preferred 
direction  like  height  in  the  atmosphere  or  depth  in  the  earth, 
a useful  approximation  is  to  regard  the  medium  as  having  a 
layered  structure.  Thus,  if  the  regular  term  <x  (?)  could  be 
neglected  in  each  layer,  we  could  simply  write  for  (2.1.1.) 


a (r)  » an  + e (?) 


where  we  now  have 


E{*(r) } = E{aQ } + E{e(?) } = a0 


(2.1.3) 


(2.1.4) 


The  random  inhomogeneities  are  thus  interpreted  as  fluctuations 
round  the  -expected  value  of  the  parameter  considered.  Follow- 
ing the  theory  of  Chernov  (I960)  for  ac-astic  waves  in  a random 

medium,  we  shall  adopt  this  simplified  case  as  a starting 
point. 
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2.2 


Description  of  a Random  Medium 

It  should  be  recalled  that  eqs.  (2.1.1)  and  (2.1.3)  are 
stationary  versions  of  time  dependent  forms.  However,  when 
changes  in  the  properties  of  the  medium  are  slow  compared  v ith 
the  duration  of  phenomena  to  be  studied,  a time  independent 
formulation  is  applicable.  In  the  following  we  shall  let  a 
denote  P wave  velocity,  but  the  random  medium  description 
by  the  use  of  a will  apply  to  the  density  p as  well  as  the 
elastic  parameters  A and  p. 


Introducing  the  refractive  index 

a. 


n (r)  = 


0 


a (r) 


and  its  deviation  from  the  mean  value  of  unity 


H(r)  = n(r)  - 1 


(2.2.1) 


(2.2.2) 


we  describe  a medium  with  random  inhomogeneities  by  considerinc 
these  functions  as  random  in  space.  The  fluctuations  in  wave 
velocity  are,  of  course,  caused  by  fluctuations  in  density 
and  the  Lam§  constants. 


Since  inji  random  medium  the  fluctuations  in  the  refractive 
index  n(r)  are  not  known,  we  shall  assume  that  its  functional 
description  represents  a random  process  (Bendat  and  Piersol 
1966).  The  variation  of  the  index  of  refraction  within  a 
volume  or  along  a particular  path  is  considered  to  be  a 
realization  or  sample  function  from  an  ensemble  of  possible 
functions  (n<r>)  ({)  denotes  ensemble).  The  realization  is 
considered  to  correspond  to  a point  in  a suitable  probability 
space  so  that  it  is  meaningful  to  speak  about  statistical 
or  ensemble  averages.  Thus  in  order  to  characterize  a random 
medium,  we  have  to  use  statistical  concepts  where  space  plays 

the  role  of  time  in  the  general  definition  of  erqodic  random 
processes . 
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We  shall  assume  that  the  random  medium  is  spatially  homogeneous 
and  statistically  isotropic.  By  spatial  homogeneity  we  mean 
that  whenever  v is  a subvolume  of  the  random  medium  and  r is 
the  fixed  difference  between  possible  point  locations  in  v 
with  position  vectors  r'  and  r'  + r then 


E{fT(?)}  = E{[n(r')r,  (?«+?)  ] } = N'(?) 


(2  7.3) 


where  the  bar  denotes  averaging  over  all  r'  within  v.  It  is 
re^ily^nderstood  that  N'  (r)  is  defined  as  the  average  product 
n(r')n(r'+r)  taken  over  the  entire  medium.  Statistical  isotropy 
eliminates  the  dependence  on  the  direction  of  r so  that 


where  r denotes  the  modulus  of  the  spatial  difference  vector 
r.  This  implies  that  when  |r|  =0 


requiring  that  the  mean  square  fluctuations  in  the  refractive 
index  be  constant  throughout  the  random  medium.  In  this  case 
we  define  a correlation  coefficient  by 

N(r)  = N'(r)/N'(0)  = N'(r)/  n2  [2.2.6) 

As  pointed  out  by  Chernov  (I960),  there  seems  to  be  no  simple 
way  to  determine  this  correlation  function  theoretically  and 
a proper  starting  point  is  therefore  some  empirical  correlation 
function.  Experience  shows  that  a large  number  of  estimated 
correlation  functions  for  fluctuations  of  physical  parameters 
in  nature  approximately  satisfy  a correlation  coefficient 
exponentially  decreasing  with  distance  like 


N (r)  = o~r/a 


(2.2.7) 


where  a is  constant.  However,  it  can  be  shown  (Chernov,  1960) 
that  the  assumption  about  n (?)  being  a continuous  function 
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implies  that 


dN  (r) 
dr 


r=0 


0 


(2.2.8) 


The  correlation  function  given  in  (2.2.7)  does  not  satisfy 
this  condition,  and  a somewhat  different  function  satisfying 
(2.2.8)  has  to  be  chosen.  A definition  chosen  by  Chernov 
(1960)  for  mathematical  convenience  which  will  be  assumed  in 
the  following  is  the  Gaussian  function 

N (r)  = e r / (2.2.9) 

where  a is  known  as  the  correlation  distance.  This  function 
is  depicted  in  Fig.  2.1.  The  correlation  distance  corresponds 
to  the  point  where  the  correlation  function  (2.2.9)  is  down 
to  e =0.37,  and  the  statistical  dependence  between  fluctuations 


Fig.  2.1  Gaussian  correlation  function  assumed  for  the  refractive 
index  fluctuations. 
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becomes  small  compared  to  unity.  Physically  the  correlation 
distance  is  interpreted  as  the  average  size  of  the  inhomo- 
geneities in  the  medium,  i.e.,  the  diameter  of  the  volume 
throughout  which  we  can  assume  the  physical  conditions  to 
be  constant. 


2 . 3 The  Larth  as  a Random  Medium 

Modelling  the  earth,  or  part  of  it,  as  a random  medium  ac- 
cording to  the  definition  in  the  previous  two  sections  is  a 
question  of  considerable  difficulty.  Due  to  the  known  varia- 
tion of  wave  velocity  with  depth  as  well  as  lateral  variation 
in  this  quantity,  only  isolated  layers  or  volumes  would  per- 
fectly pertain  to  the  random  medium  concepts  homogeneity 
and  isotropy.  Capon  (1974)  suggested  to  solve  this  problem 
by  modelling  layered  structure  of  the  earth's  crust  and  upper 
mantle,  with  P-velocity  fluctuating  around  a constant  value 
within  each  layer.  He  also  discusses  the  possibility  of  de- 
comp Dsing  this  model  into  an  equivalent  deterministic  and 
an  ideal  random  medium,  where  the  random  medium  only  gives 
rise  to  fluctuations  in  amplitude  and  phase.  Such  a model  is 
realistic  when  the  angles  of  incidence  at  the  interface";  are 
small  for  all  portions  of  the  wave. 


The  earth  model  assumed  throughout  this  text  is  the  one 
depicted  in  Fig.  2.2,  where  random  inhomoqeneities  are 
confined  to  the  upper  part  of  the  mantle  and  the  crust.  It 
is  obvious  from  this  figure  that  the  'histories'  of  the  two 
rays  in  the  inhomogeneous  medium  on  the  source  side  of  the 
path  are  about  identical,  ensuring  a plane  wave  to  enter  into 
the  inhomogeneous  medium  on  the  receiver  side.  This  assumption 
is  important  in  the  application  of  the  theory  to  be  developed. 


r 


Fig.  2.2  Schematic  earth  model  with  inhomogeneities  confined  to  the 

upper  50C  km  (beyond  dotted  line) . The  rays  drawn  correspond 
to  the  extreme  ones  for  NORSAR  (6Awl  deg)  . 
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3-  SCATTERING  Or  P- WAVES  IN  A RANDOM  MEDIUM 
^ ^ The  Equation — of  Motion  In  a Weakly  Inhomogeneous  Medium 

The  only  means  of  solution  of  the  linear  elasticity  relations 
in  an  inhomogeneous  medium  is  through  approximation  methods 
which  impose  restrictions  on  the  structural  complexity  of 
the  medium.  However,  even  crude  approximations  are  justified 
through  the  practical  usefulness  of  exact  of  asymptotic 
solutions  of  the  wave  equation.  Assuming  infinitesimal  strain 
in  a perfectly  elastic  and  isotropic  solid  where  u(r,t)  denotes 
displacement  about  the  equilibrium  position  of  a point  P(r) 
we  obtain  from  Karal  and  Keller  (1964),  applying  the  operator 
equality  V2u  = V (V*u) -VxVxu: 

32u 

P ~T  (A+2u)  V(V*u)+VAV*u-ljVxVxu+Vux(7xu)+2  (Vp*V)u  (3.1.1) 

where  Mr),  U (r ) and  p(r),  the  elastic  moduli  and  density 
respectively,  are  random  functions  in  space  according  to 
eq.  (2.1.3).  Theoretical  seismology  usually  deals  with  wave 
propagation  through  media  where  gradients  of  the  medium's 
characteristics  are  either  smooth  and  small  or  abrupt,  in 
both  cases  solution  of  eq.  (3.1.1.)  involves  two  distinct 
waves,  dilatational  and  shear.  The  amount  of  coupling  between 

these  waves  depends  on  the  above-mentioned  gradients  (Landers, 
1971) . 


The  case  to  be  treated  here  is  a weakly  inhomogeneous  medium, 
where  the  spatial  change  of  the  elastic  parameters  is 

small  as  compared  to  the  parameters  themselves.  In  this  case 
we  simply  neglect  terms  containing  gradients  of  the  medium's 
characteristics.  Rewriting  eq.  (3.1.1)  using  the  expression 


for  dilatation,  we  obtain  the  modified  form 
92u 

P =( A + 2u ) VO-pVx ( Vxu) 

at2 
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Taking  the  divergence  of  this  equation,  remembering  that  terms 
containing  Vp,  VX  and  Vp  are  negligible  gives: 


Dividing  through  by  X + 2p  we  arrive  at  the  well-known  differential 
equation: 


-1  ill 

a1 2  9t2 


v2  e 


o 


(3.1.5) 


where 


a = a (r)  = 


X (r) + 2p (r) 
P (r) 


(3.1.6) 


is  the  velocity  of  the  primary  disturbance  expressed  by  0. 

If  the  amount  of  variation  in  the  elastic  parameters  is  small, 
say  a few  per  cent  about  the  mean  value,  this  property  will 
certainly  also  be  shared  by  the  wave  velocity  a (cf.  eq.  (3.1.6)). 
Introducing  n as  defined  in  eq.  (2.2.2)  we  infer  that 
(see  Chapter  2) : 

1 n2 


where 

(3.1.8) 

The  wave  equation  for  a weakly  inhomogeneous  medium  then  finally 
becomes : 


_ ( 1+n) 


a 


(3.1.7) 


(1+n)2  9 2 0 


a 


9t‘ 


- V2  0 = 0 


(3.1.9) 
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This  is  the  form  of  the  wave  equation  for  compressional  elastic 
waves  required  in  order  to  take  advantage  of  the  simplicity 
of  Chernov  s theory  for  acoustic  waves.  Note  that  since  n is 
an  unknown  function  of  position,  there  is  no  easy  way  to  solve 
the  differential  equation  (3.1.9). 


3 . 2 Rytov's  Method 

Random  inhomogeneities  cause  fluctuations  both  in  amplitude 
and  phase.  This  is  a statement  based  on  practical  experience 
(fading  of  radio  signals,  twinkling  of  stars)  as  well  as 
theoretical  considerations.  Imagine  a uniform  plane  mono- 
chromatic wave  of  P-type  in  homogeneous  medium  travelling 
into  a weakly  inhomogeneous  medium  of  the  type  described  in 
2. 1-2. 2 (see  Fig.  3.1).  The  primary  waves  in  the  inhomogeneous 
medium  are  assumed  to  satisfy  eq.  (3.1.9).  In  the  homogeneous 
medium  we  write  the  plane  wave  in  the  standard  form: 


where  for  simplicity  the  direction  of  the  x-axis  coincides 
with  the  normal  to  the  wavefront. 

Due  to  the  effect  of  the  inhomogeneities,  eq.  (3.2.1)  has  to 
be  modified  accordingly.  The  amplitude  and  spatial  part  of  th£ 
phase  Lunction  are  now  supposed  to  be  unknown  functions  of 
position  in  the  inhomogeneous  medium,  and  we  are  looking  for 
a solution  of  eq.  (3.1.9)  of  the  form: 

0 = A (r)  e"i(u)t_S{^)  } 

' ’ e (3.2.2) 


The  amplitude  function  A(r)  is  alternatively  expressed 
as 

A (? ) = Aq  e l°g  (A  (r  )/Aq) 

where  log  means  the  natural  logarithm  to  be  taken. 


Fig.  3.1  Schematic  illustration  of  forward  scattering  in  an 

inhomogeneous  medium.  The  scattering  field  at  a point 
P(r)  is  given  as  the  integrated  effect  within  a cone 
with  vertex  in  P and  aperture  angle  1/ka  {the  Fresnel 
approximation) . 


Inserting  this  expression  into  eq.  (3.2.2)  we  get 

-iUt-S(r)+i  log  -i  (ut-Y  (r)  ) 


9 = A e 
o 


= Ao  e 


(3.2.4) 


where 


'J'(r)  = S (r)  - i log  (Mil) 

A 

o 


(3.2.5) 


The  real  part  of  the  complex  function  Y is  the  spatial  phase 
function  S (r ) , and  the  imaginary  part  expresses  the  natural 


logarithm  of  the  ratio  between  amplitudes,  log  (A(r)/A  ).  Note 


that  the  latter  ratio  also  denotes  the  fluctuations  in  log- 
arithmic amplitude  (logamplitude)  relative  to  the  incident 
wave.  The  idea  now  is  to  use  the  right  hand  part  of  eq.  (3.2.4 


iii  the  basic  wave  equation  for  the  inhomogeneous  medium 
(eq.  (3.1.9)).  Doing  this,  we  obtain  an  equation  where 
0 has  been  replaced  by  Y,  namely 


(W)  2-iV2f  = (1+n/2 

a 2 
o 

In  the  homogeneous  medium  the  function  4*  has  the  form 

*o=kx  according  to  its  definition,  and  satisfies  the 
equation 


(V4'q)  2-iV2t'o  = k2 


(3.2.7) 


where 


is  the  wavenumber  of  the  homogeneous  wave. 


(3.2.8) 


Writing 


V = 4*  +4', 

o 1 


(3.2.9) 


we  have 


’•'l  ~ si  ~ i log  (a/ao) 


(3.2.10) 


and 


S1  = S " kX  (3.2.11) 

We  are  now  in  a position  to  successfully  combine  eq.  (3.2.6) 

and  eq.  (3.2.7).  Invoking  eq.  (3.2.9)  into  (3.2.6)  and  sub- 
tracting eq.  (3.2.7)  yields: 


m"  9 ^ 


r* 


l I 'III  | J P^fVM 
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2W0V1'1  - iv2^  = 2nk2  +[ n2k2-  (v^) 2 ] 


(3.2.12) 


The  next  step  is  to  drop  the  terms  in  the  square  brackets 
under  the  assumption  about  small  perturbations  in  the  re- 
frac-ive  index,  so  that  nz  can  be  neglected  compared  to  n. 
The  additional  requirement  is  that  and  according  to 

(3.1.8)  this  imposes  the  restriction  that 

El7,t,ll  « 1 (3. 


Eq.  (3.2.12)  is  then  simply  written  as  the  linear  differential 
equation . 


2V4'q74'1  -iV2  4'1  = 2nk2 

Remembering  that  4*  =kx,  we  have 

o 

2k(34'1/3x)  - iV2t1  = 2n k 2 
and  by  defining  a new  function  W related  to  4^  as 
W = 4'1  eikx 

we  obtain  the  inhomogeneous  Helmholz  equation 
V2w  + k2w  = i2nk2elkx 
for  the  function  W. 


(3.2.14) 


(3.2.15) 


(3.2.16) 


(3.2.17) 


Eq.  (3.2.14)  was  obtained  under  the  restriction  that  the 
phase  change  and  relative  change  in  amplitude  per  wave  length 
is  small  (cf.  eq.  (3.2.13)).  Because  this  essentially  is  a 
limitation  on  the  gradient  of  4.,,  it  obviously  serves  as  a 
better  approximation  than  the  method  of  small  perturbations 
which  assumes  that  the  perturbed  quantity  is  small  (Born 
approximation).  However,  neglecting  (V4^)2  while  retaining 
V24^  in  eq.  (3.2.12)  is  a questionable  step  which  has  been 
commented  on  by  Taylor  (1967).  As  far  as  can  be  seen,  there 
is  no  a priori  reason  why  the  two  terms  in  the  square  brackets 
of  eq.  (3.2.12)  simply  cancel  each  other  and  thus  justify 
an  assumption  of  the  type  above. 

The  applicability  of  Rytov's  method  turns  out  to  be  conditioned 
upon  the  inequalities: 

£|V  log (^—) I <<  1 (3.3.1) 

and 

<<  1 (3.3.2) 

Inequality  (3.3.1),  which  states  that  the  scattered  energy 
in  going  a wave  length  be  small,  is  plausible  in  view  of  the 
restriction  |n|<<l.  Since  we  write 

S = SQ  + Sl  = kx  + S:(r)  (3.3.3) 

the  inequality  (3.3.2)  may  alternatively  be  written  as 

£|VS|  ~ 1 (3.3.4) 

Expressing  the  gradient  of  S as  components  in  a rectangular 
coordinate  system  (see  Fig.  3.1),  we  have 
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I 3S|  L. 

'W  ~ k' 


9S 

ay 


V o 

377  <<  k, 


as, 

az  i 


<<  k 


(3.3.5) 


These  expressions  show  that  the  transverse  components  of  the 
gradient  of  the  phase  (i.e.,  the  transverse  wave  wavenumber 
components,  are  small  compared  to  the  longitudinal  component. 

is  means  that  the  direction  of  propagation  for  waves  in 
the  inhomogeneous  medium  deviates  only  slightly  from  the 
original  direction  in  the  homogeneous  medium. 

Sharply  directed  scattering  is  known  to  be  produced  by  large 
scale  inhomogeneities  where  the  condition 


ka  >>  1 


(3.3.6) 


IS  fulfilled  (Morse  and  Feshbach,  1953).  In  inhomogeneous 
media  of  this  type  we  can  neglect  back-scattering  and  wave 
ef lection  in  the  sense  that  the  energy  lobe  for  a point 
source  on  the  wave  front  will  have  the  form  illustrated  in 
Fig.  3.1.  This  serves  as  a simplification  for  the  solution 
OF  eq.  (3.2.17)  at  a point  P(?)  in  the  inhomogeneous  medium. 

e actual  solution  of  eq.  (3.2.17)  is  a Green’s  function 
problem  (Morse  and  Feshbach,  1953)  and  the  general  form  is- 


w = - 


4 TT 


f -i2r^k2  pi  (k  u+x’ ) ) 

l e dv 


(3.3.7) 


where jc  r-r  | in  accordance  with  the  symbols  used  in  Fig.  3, 

is  e part  of  the  inhomogeneous  medium  contributina  siani- 

Urge  T S°1Uti°n  **  P°int  P‘  The  sumptions  about 
rge  scale  inhomogeneities  justify  the  limitation  of  the 

vo  ume  V to  the  layer  which  lies  between  the  base  of  the 
inhomogeneous  medium  (x-0  in  Fig.  3.1)  and  P,  and  within  a 
con©  with  its  VGirtGx  p ar^A  a 

, n/.  , P and  an  aPerture  angle  of  the  order 

1.0/ka  (see  Fig.  3.1) . 


Within  these  limits  we  approximate  the  length 


i = ( (x-x* ) 2 + d2 ) ^ 


(3.3.8) 
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8.  ~ x-x'  + i 2 , 
2 x-x ' 


(3.3.9) 


where 


d2  = (y-y ' ) 2 + (z_z i ) 2 


(3.3.10) 


Multiplying  eq.  (3.3.7)  by  e~lkx  we  obtain  the  solution  for 
the  function  Y yielding: 


f - - lk  r n Q 

1 2 7T  ^ Z e 


n ik ( £- (X-x1  ) ) 


(3.3.11) 


Replacing  1.0/i  by  1.0/x-x'  and  substituting  in  the  exponent: 
of  eq.  (3.3.11)  from  eq.  (3.3.9)  gives 


ii/  - _ tk2  , n 

1 2 IT  l e 


-d  } 
( x— x * ) ’ 


(3.3.12) 


Taking  real  and  imaginary  parts,  we  finally  obtain  the  ex- 
pressions for  the  scattered  field: 

S.  = &L  f f T -sin(kdz/2  ( x— x * ) ) , , 

1 2tt  i 1 I x^ n (x1  ,y ' , z ' ) dx  'dy  'dz  ' (3.3.13a! 


log  A/A  = / / 7 cos (kd2/2 (x-x1 ) ) , 

° 2tt  Jq  l l h (x ' ,y ' , z ) 

dx ' dy 1 dz ' 

* (3.3. 13t 

These  equations,  which  turn  out  to  correspond  to  the  Fresnel 
approximation  in  diffraction  theory  (Jenkins  and  white,  1957), 
are  the  starting  point  for  the  statistical  derivations. 
Noteworthy,  eqs.  (3.3.13)  may  be  thought  of  as  the  analytical 
expressions  foi  the  effects  causing  crinkled  wavefronts  to 
P-waves.  In  the  next  chapter  these  effects  will  be  implemented 
in  a refined  model  for  travel  time  and  logamplitude  observations. 


"^TT~ — 
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3.4  Statistical  Quantities  of  the  Scattered  Wave  Field 

Since  a detailed  derivation  of  the  statistical  aspects  of 
the  wave  scattering  problem  has  been  carried  out  by  Chernov 
(I960),  only  main  formulae  will  be  summarized  and  important 
steps  stressed  in  this  section. 

Assuming  that  the  receiver,  i.e.,  our  seismograph,  is  located 
at  the  point  Q'  with  spatial  cartesian  coordinates  (L,0,0) 

(cf . Fig.  3.1),  our  first  aim  is  to  obtain  expressions  to 
the  mean  square  fluctuations  in  phase  and  logamplitude . 

Using  the  symbol  4>  for  the  phase  fluctuations  S.,  and  A for 
the  logamplitude  fluctuations  log (A/Aq) , we  rewrite  eqs . 

( 3 . 3 . 13a ,b)  as : 

L + 00 

$ = / / / F . (L-x ' , k , d ) g (x ' ,y' , z ' ) dx ' dy 1 dz ' (3.4.1a) 

0 - * 

L + oo 

A = / / / F_ (L-x' ,k,d) n (x' ,y ' ,z' )dx'dy 'dz ' (3.4.1b' 

0 - 00  C 

where  F ^ and  F2  are  easily  deduced  from  eq.  (3. 3. 13a, b). 
Squaring  each  of  these  expressions  and  taking  the  average, 
remembering  from  section  2.2  that 

n(r')n(r")  = N'(|r'-r"|)  = n2N(r)  (3.4.2) 


where 


r2  = (x'-x")2  + (y'-y")2  + (z'-z'')2  (3.4.3) 

we  have 

— — LL  + 00 

$2  = n2  II  llll  F (L-x' ,k,d’)F.  (L-x"  ,k,d")N(r)  (3.4.4a) 

00  - * 1 1 


dx ' dy ' dz ' dx ' ' dy ' ' dz ' ' 
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A2 


— LL  + oo 

= n2  //  ////  F ~ (L-x ' ,k ,d 

00  - oo  2 


)F2 (L-x' • ,k,d* ' )N(r) 


(3.4.4b) 


dx ' dy ' dz ' dx ' ' dy • ' dz • • 


Using  linear  coordinate  transformations  and  variable  shifts, 
we  obtain  under  the  assumption  of 


that 


F = k-aMl 

2 

A 2 * n 2 k2  a L ( 1 


arctan  D) 
arctan  D) 


provided 


(3.4.6a) 

(3.4.6b) 


(3.4.7) 


a dimensionless  quantity  denoted  wave  parameter,  is  inter 
mediate  or  large  compared  to  unity.  Note  that  the  case 
D<<1  is  excluded.  An  important  result  is  that  the  ratio 
RMS  (root  mean  square)  logamplitude  to  RMS  phase  given  by 


aA  .A2 . ^ 


1-  p arctan  D 


L 1+  q arctan  D 

is  a quantity  always  less  than  unity. 


(3.4.8) 


A more  complete  characterization  of  the  statistical  properties 
of  the  wave  field  is  achieved  by  using  correlation  ^unctions. 
The  cross  correlation  between  logamplitude  and  phase  fluctua- 
tions analogous  to  (3.4.4)  is  given  by: 
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— LL  + 


^ n 00  Fl(L"X,'k'd')F2(L-x,'^,d")N(r 

dx'dy'dz 'dx' ’dy ' 'dz ' ' 


) (3.4.9) 


Applying  linear  transformations  and  variable  shifts  and 
assuming  as  before  L,>a,  eg  (3  4 <»  „ . ..  r d 

or  large  compared  to  on^  ^ **  ° 


TK  = 


/tT  ~ , 


rg-  k2  a3  log  (1+d 2 ) 


(3.4.10 


We  are  now  able  to  find  an  expression  for  the  correlation 

atVpointV (LtoeO)°9arithmiC  a,"PlitUde  and  Phase  fluctuation! 
point  Q a, 0,0).  usr„g  for  this  quantity  the  definition: 


' <J>A 


4>A 


y/TT  /7C' 


(3.4.11 


we  have 


‘♦A 


= 1 log  ( l+p ) 2 

^u1- (arc tan  D)  1 


(3.4.12 


!rthethaV£  “e  COnSider  ° aS  3 raeaSUre  £°r  the  linear  extent 
at  sJlind°:rne°US  medi“'  the  correlation  which  exists 
10,  :;“  VanUheS  “ ° — .M  approaches 

The  next  step  is  to  consider  the  question  of  autocorrelation 

;s:r;hr  io9amputude)  ■*  «««•■*  ^in,  po^s 

Assume  that  our  receivers  are  located  along  the  direction 

a :r?r Let  the  c°°rdinates  °£  «»  -mv-  b; 

' (L+iL'°'0)  corresPonding  to  o'  and  Q"  ln  Fig  3 l 

::  ;;rral  ^ 
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VAL>  = 


L 

/ 

0 


L+AL  + <» 

/ mi 

0 — 00 


F1 (L-x'/k,d,)F1 (L+AL-x ' • , k , d ' ' ) 


N (r)dx'dx ' *dy 'dy ' 'dz'dz' ' (3.4.13 

It  should  be  noted  that  because  of  the  directional  character 
of  the  scattered  energy,  waves  scattered  in  the  layer  bounded 
by  the  planes  x=L  and  x=L+AL  aie  incident  on  the  second  but 
not  on  the  first  receiver.  Therefore  these  waves  can  be  neq- 
lected  in  calculating  the  autocorrelation  replacing  the  limits 
corresponding  to  the  integration  variable  dx ' by  L instead  of 
L+AL.  Variable  shifts,  linear  tranf ormations  of  origin  and  the 
condition  L>>a  justify  the  formula 


r^AL)  = Ra(AL) 


1 + ( 


?AL}  2 
ka2 


(3.4.14) 


where  ( AL)  is  the  longitudinal  autocorrelation  for  log- 
amplitude  fluctuations  defined  as  eq.  (3.4.13)  with  F 

replaced  by  ?2 . Eq.  (3.4.14)  was  obtained  under  the  assumption 
that  D is  large  compared  to  unity. 

For  large  values  of  D we  can  write  eqs.  (3.4.6)  in  the  form 


72  = J2 


L 


(3.4.15) 


Defining  the  correlation  coefficients 


r<*>  ( AL) 


VAL) 


we  have 


Ra (AL) 
r.(AL)  = 

I2 


(3.4.16) 


r$(AL)  = rA(AL)  = ± 

1 + (lAL) 2 

ka2 


(3.4.17) 


This  formula  shews  that  if  Al^a,  then  the  fluctuations  in 
the  longitudinal  direction  are  almost  completely  correlated 
in  a large  scale  weakly  inhomogeneous  medium  (ka>^l).  De- 
fining the  longitudinal  correlation  distance  as  the  distance 
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ALr  where  the  correlation  coefficient  is  down  to  0.5,  we  have 


al.  = Sai 

k 2 


(3.4.18) 


Eq-  (3*4*17)  corresponds  to  a correlation  function  quite 
commonly  used  for  representing  the  anomalous  gravity  field 
of  the  earth  (Heiskanen  and  Moritz,  1967). 

Finally  we  look  at  the  expressions  for  the  transverse  auto- 
correlation of  the  amplitude  and  phase  fluctuations.  We  shall 
let  the  receivers  be  located  tangentially  to  the  wave  front 
like  P and  Q'  with  coordinates  (L,y,0)  and  (L,0,0)  respectively. 
Expressions  for  the  fluctuations  are  still  given  by  eqs . 

(3.4.1)  replacing  d by  d ' = /(F-y'r+z'2  for  receiver  P and 
^ = •V  ' * +z ' ' z for  receiver  Q' . 

The  transverse  autocorrelation  for  phase  fluctuations  is 
given  by: 

L L + oo 

Vy)  = n?0  0 Fl(L"X''k'd,)Fl(L-x’'k^")N(r)  (3.4.19) 

dx'dx' 'dy'dy ' 'dz 'dz ' ' 

RA(y),  the  transverse  autocorrelation  for  logamplitude , is 

analogously  given  by  eq.  (3.4.19)  replacing  F,  by  F 

1 J 2 ’ 

Defining  normalized  functions  as  before 

r fvi  - R$(y)  (y) 

r<*,(y)  = - - A 


r.(y)  = 


(3.4.20) 


Chernov  has  shown  to  obtain  for  L>> 
the  formulae: 


a and  large  values  of  D 


(3.4.21a) 


r*(y>  = 


'y/a  -5<7-si'5^»> 


1 - - arctan  D 


rA(y)  - 


*~y/a 

i + I 

D arctan  D 


where 


Si (x)  = 


/ sis_t  dt 


(3.4 


(3.4 


The  graphs  of  the  functions  (3.4  21a, b)  with  D=10  are  shown 
in  Fig.  3.2  compared  to  the  correlation  functions  for  the 
fluctuations  in  the  refractive  index.  Note  that  because  of 


Fig.  3.2  Transverse  autocorrelation  function  for  phase  (3)  and 

logamplitude  (1)  together  with  the  Gaussian  correlation 
function  (2).  Redrawn  from  Chernov  (1960). 


21b) 


22) 
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the  rotational  symmetry  assumed,  y should  be  interpreted 

as  the  distance  r between  receivers  tangentially  to  the 

wave  front  ra.her  than  the  coordinate  in  a specified  direction. 

The  important  result  to  stress  is  that  the  transverse  auto- 
correlation functions  for  logamplitude  and  phase  extend  to 
a distance  of  the  order  of  the  correlation  distance  in  the 
medium.  Recently,  Christof fersson  (1974)  has  pointed  out 
that  eqs . (3.4.21)  are  valid  for  y/a<  1 since  Chernov  did  the 
approximation 


e ~ 1 (3.4.23) 

to  obtain  these  equations.  Christof fersson  (1975)  suggests 
using  a series  expansion  in  the  actual  approximations  to 
obtain  transverse  autocorrelation  functions  valid  for  a wide 
range  of  y/a  and  D values.  In  this  work  no  attempt  has  been 
made  to  correct  eqs.  (3.4.21)  for  y/a  values  considerably 
larger  than  1,  and  this  should  be  remembered  below. 


3 . 5 Classification  by  the  Wave  Parameter  D 

Seismological  wave  propagation  problems  fall  into  three 
main  categories  depending  on  the  wave  parameter 


Limiting  the  discussion  to  the  case  ka>>l,  i.e.,  large  scale 
inhomogeneities,  we  have: 

1)  The  region  D<<1  where  the  Uniform  Plane  Wave  Model  can 
be  used.  As  the  extent  L increases  and  approaches  the 
correlation  distance  in  size,  we  have  to  account  for 
the  variation  in  the  physical  properties. 
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2)  When  Del,  Ray-Theory  might  be  applicable. 

3)  The  region  D>1  where  Diffraction-Theory  is  necessary 
and  the  stochastic  behavior  of  the  wave  field  should 
be  regarded  as  significant. 

The  preceding  theory  treats  P and  S as  uncoupled  waves,  where 
P energy  is  converted  into  P energy  due  to  the  inhomogeneities 
in  the  medium.  Knopoff  and  Hudson  (1964,  1967)  show  that 
P-S  and  S-P  conversion  can  safely  be  neglected  compared  to 
P-P  under  circumstances  when  ka>>l.  The  important  thing  to 
remember  is  that  deviations  from  the  weakly  inhomogeneous  case 
certainly  exist  in  the  real  earth,  and  that  these  deviations 
may  lead  to  considerable  P+S  conversion  as  well  as  back- 
scattering  effects. 

Point  (3)  in  the  previous  classification  scheme  requires 
a statistical  attitude  to  seismic  problems.  An  important 
assertion  put  forward  and  made  use  of  in  the  following  is 
the  normality  in  the  distribution  of  the  random  variables 
* and  A.  In  fact,  if  the  entire  distance  L which  contributes 
to  the  values  of  4>  and  A is  divided  into  segments  of  the 
order  of  the  correlation  distance,  then  the  changes  in  <t> 
and  A along  different  segments  will  be  almost  statistically 
independent.  According  to  the  Central  Limit  Theorem  of 
probability  theory,  when  the  number  of  independent  segments 

is  large,  quantities  like  4>  and  A will  obey  the  normal  dis- 
tribution law. 
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4.  THE  APPLICABILITY  OF  CHERNOV ' S MODEL 

__ 

4 . 1 Selection  of  Data 

Normally,  frequency  domain  analysis  of  data  is  preferred  when 
formulae  obtained  involve  the  wavenumber  k as  in  Chapter  3. 

On  the  other  hand,  using  spectra  calculated  from  a data  window 
according  to  some  estimated  or  assumed  slowness  introduces 
at  least  two  error  sources.  First  of  all,  back-scattering 
or  coda  effects  at  the  end  of  the  window  is  likely  to  play 
a role.  Secondly,  the  estimate  of  window  start  time  is  asso- 
ciated with  some  bias  or  measurement  error.  Time  domain  analysis 
was  therefore  regarded  as  favorable  in  this  case.  Narrow 
bandpass  filtering  was  obtained  by  exposing  the  20  Hz  sampled 
signals  to  a Butterworth  3rd  order  lowpass  filter  followed 
by  a similar  highpass  filter.  The  combined  filter  shows 
a very  distinct  peak  response  at  the  average  frequency  and 
was  advantageous  compared  to  conventional  narrow  bandpass 
filters  with  a reasonable  number  of  bandpass  coefficients. 

For  analysis  purposes  10  teleseismic  P events  were  chosen 
as  the  basic  data  source.  Other  selection  criteria  were 
good  signal-to-noise  ratio  as  well  as  coverage  in  azimuth 
and  velocity.  NOAA  epicenter  information  for  all  events 
is  given  in  Table  1. 

It  was  decided  to  measure  travel  time  and  amplitude  as  soon 
as  possible  after  the  onset  of  each  event,  thereby  eliminating 
as  much  as  possible  of  the  back-scattering  and  coupling  effects. 
From  a computer-plotted  display  of  the  filtered  traces  the 
peak  (or  trough)  of  the  first  (or  second)  cycle  was  uniformly 
identified  by  visual  inspection.  At  the  same  time  sensors 
with  bad  performance  due  to  communication  problems  were  masked 
in  order  to  avoid  abnormal  data  in  the  computer  programs. 

The  exact  travel  time  and  amplitude  for  the  identified  phase 
were  determined  as  the  extremum  of  a parabola  fitted  to  seven 
points  around  the  peak.  The  accuracies  in  these  time  and 
amplitude  determinations  are  not  easily  fiaured  out  since 
the  exact  value  of  what  we  are  trying  to  measure  is  unknown. 


NOAA  EPICENTER  INFORMATION 


Region 


Event 

No. 

1 

(*> 

2 

(*) 

(***) 

3 

(*) 

4 

(*) 

(***) 

5 

<*) 

<**) 

6 

(*> 

- 

C***) 

7 

(*) 

8 

(*) 

9 

<*> 

10 

(*) 

(**) 

Date 

(m/d/y) 

11/21/72 

12/10/72 

OS/ 17/ 7 3 
00/06/72 

04/10/72 

04/11/72 

07/27/71 

07/00/72 

07/03/71 

03/26/71 


Origin  Time 
(h  n s) 

17.01.55. 3 
03. 35.48.2 

09.30.00. 9 

01.12.50.4 

15.07.49. 1 
02.21. 15.7 

02.02.49.6 

12.10.54.7 

14.00. 00. 1 
17.35.18.0 


Aleutians 

Japan 

Kazakh 

Iran 

Tanganyika 

Mid-Atlantic 

Peru 

Mexico 

Nevada 

Alaska 


TABLE  1 

Events  used  in  analysis 

(*)  Events  filtered  at  0.7  Hz  used  in  estimating  Chernov  (1960) 
parameters  and  in  slowness  estimation  problems 

(***)  Events  filtered  also  at  0.9,  1.1  and  1.3  Hz. 

See  references  to  these  symbols  in  the  text. 
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Fig.  4.1  Determination  of  exact  travel  time  and  amplitude  (open 

circle)  by  fitting  a parabola  to  seven  points  around  the 
pea  sample  (black  dots).  Sampling  interval  is  0.05  sec. 


A few  examples  of  determination  of  peaks  by  this  method  are 
shown  in  Fig.  4.1,  from  which  we  conclude  that  the  measured 
times  are  certain  within  a fraction  of  a sampling  interval 
of  50  milliseconds.  Of  course,  this  conclusion  is  based 
on  an  assumption  about  sufficient  coherency  in  the  signals. 
All  the  observed  times  were  reduced  to  common  observation 
plane  by  correcting  for  differences  in  altitude  among  the 
seismic  stations.  This  correction  was  done  using  the  ap- 
parent velocity  determined  by  the  NORSAR  Event  Processor 
(EP)  G j^ystdal  (1973).  Finally,  it  should  be  mentioned  that 
the  EP  determined  velocity  and  azimuth  were  always  used 
when  projecting  horizontal  distances  between  sensors  into 
the  wave  front  to  obtain  the  transversal  distance  (PQ' 
as  opposed  to  PQ' • in  Fig.  3.1). 
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4.2 


Observations 


The  validity  of  the  principles  outlined  in  this  chapter  will 
be  quite  general,  although  the  headline  tries  to  relate  the 
content  to  a special  topic,  namely,  seismic  measurements. 

Since  the  data  analysis  in  the  following  will  be  carried  out 
for  a network  or  an  array  of  stations,  there  will  be  a need 
for  multidimensional  concepts  in  terms  of  vectors  and  matrices. 
In  this  specific  context  T will  denote  a column  vector  of 
n observations.  The  quantity  expressed  by  T will  be  either 
phase  travel  times  or  logamplitudes  as  observed  by  n dif- 
ferent sensors  for  short  periodic  compressional  waves.  The 
expected  value  of  the  vector  T is  the  vector  of  expectations 

of  its  elements  (Morrison,  1967).  if  E denotes  expectation, 
we  have 


E(tx} 


E { T } = 


(4.2.1) 


E { t } 
n 


where  are  the  individual  elements  of  T. 


The  concept  expectation  is  extended  epresent  any  linear 
model  formulation  in  an  appropriate  coordinate  system.  This 
means  that  we  may  write 


E{T)  = RU 


(4.2.2) 


where  R is  an  n x p matrix  of  coordinates  and  U is  a p x 1 
vector  of  p appropriate  regression  coefficients.  For  simpli- 
city we  do  not  distinguish  explicitly  between  matrices  and 
vectors . 


An  illustration  of  the  applicability  is  given  in  Fig.  4.2,  where 
the  observed  travel  time  (or  logamplitude)  field  across  an 
array  of  seismic  stations  is  schematically  decomposed  into 
a linear  deterministic  trend  effect  superposed  by  a spatially 
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STOCHASTIC  EFFECT 


^trend 
ARRAY  OIAMETER  *J 


TREND  EFFECT 


DISTANCE  ALONG  A CROSS -SECTION 


Fig.  4.2  Schematic  illustration  of  a linear  deterministic  trend 
model  superposed  by  a stochastic  scattering  effect. 


stochastic  effect  S assumed  to  satisfy: 


E{  S } = 0 


S may  be  thought  of  as  i vector  of  n elements  analog  to 
T in  eq.  (4.1.1). 


The  idea  here  is  that  the  effects  depicted  in  Fig.  4.2  are 
the  dominating  components  in  the  observed  compressional 
wave fie Id.  S corresponds  to  the  stochastic  fluctuations  given 
by  eqs.  (3.3.13),  and  which  are  physically  interpreted  as 
the  scattering  response  of  an  inhomogeneous  crust  and  upper 
mantle.  Note  that  no  measurement  error  or  noise  effect  has 
been  incluued  in  Fig.  4.1,  as  these  terms  are  thought  of  as 
oniy  producing  small  ripples  on  the  solid  line  representing 
the  total  field. 


■I 
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If  the  location  of  seismic  sensors  is  expressible  in  a local 
two-dimensional  cartesian  coordinate  system  with  x and  y 
axes  pointing  W-*E  and  S+N  respectively,  then  the  deterministic 
effect  in  the  travel  time  field  is  readily  expressed  in  terms 
of  slowness-components  (UX,U  ) (Gj0ystdal,  1973).  An 

improved  model  for  the  observed  times  is  then  given  according 
to : 


' 1 

t 

o 

t . = 
1 

rix 

• 

U 

X 

riy 

U 

y 

i- 

+ S . + £ . 

1 1 


i=l , 2 , . . . ,n  (4.2.4) 


where  (rix,  r are  the  position  vector  components  for  the 
i-th  station,  s^  is  random  scattering  effect  and  e.  is 
measurement  error  or  noise.  According  to  eq.  (4.2.1)  we  have 
from  eq.  (4.2.2)  that 


E{t.}  = E{ t +r . U +r . U +s.+e  } 
1 o IX  x ly  y i L ; 


= t +r,  U +r  U 
o ix  x ly  y 


i 1 , . . . , n 


(4.2.5) 


since 


When  t.  represents  logamplitude,  Ux  and  Uy  are  logampHtude 
gradients  in  the  x,y-plane.  No  such  gradients  are  observed 
in  real  data  and  we  may  set 


This  is  also  reasonable  in  view  of  magnitude  correction  tables 
which  allow  for  maximum  5 per  cent  variation  in  amplitude 
within  a 1 degree  teleseismic  distance  interval,  corresponding 
to  the  diameter  of  NORSAR.  For  logampl itudes  we  then  write 
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the  model 


log  A.  = log  A + a . + e 
i o i ei 


i=l. 


»n 


(4.2.8a) 


where  a.  is  the  scattering  terms  and  e measurement  error 
for  logamplitude.  Thus  we  have 


E{ log  A } = log  A 
1 ^ o 


(4.2.8b) 


If  we  let  t = loa  a a - , 

o aq,  s±  - ait  ei  = eit  we 

(4.2.4)  may  represent  time  and  logamplitude 

of  course,  with  obvious  modifications. 


see  that  eq. 
s imul taneous ly , 


In  matrix  formulation  eq.  (4.2.4)  is  written 
T = TU  + S + € 


where 


f 1 r • • • t n 


(4.2.9) 


At  this  stage  it 
notation  for  the 


is  natural  to  introduce  the  extended 
n component  vector  T,  defined  as 


variance 
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E{(T-E{T}) (T-E{T})*>  = E{ (T-RU) (T-RU) *}  = 


E{SS*}  = (4.2.10) 


'll 

J21 


nl 


2 

12 

2 

22 


2 

In 

2 


nn 


= I 


°ij  = E{<t1-E{t.})(t1-E{ti)>*}  , 1>j=1 


(star)  denotes  transponation  and  £ is  for 
included  in  the  S term. 


convenience 


The  normalized  version  of 

each  element  a2.  . by  o o 

x3  * ii  jj 


eq.  (4.2.10)  is  obtained  dividing 
yielding  the  correlation  matrix 


(4.2.11) 


The  matrix  C is  invariant 
of  the  variates  in  T. 


1 


under  changes  of  scale  and 


origin 


“erlnLtn  array  °f  -l-tivelv  closely  space  seismometers 
the  earth  s surface  to  be  located  approximately  trans- 
versa  ly  to  the  wave  front,  then  the  correlation  coefficient 
Elj  eq.  (4.1.11)  is  given  by  eqs.  (3.4.21)  orovlded 
chernov.s  ,1»60|  theory  for  wave  scatterinq  is' considered 
to  be  acceptable.  The  existence  of  matrices  Z and  C sig- 
ni  cantly  different  from  a diagonal  and  identity  matrix 
respectively  means  that  there  is  in  general  an  all-to-all 
dependence  between  our  observations.  More  precisely  formulated- 

its  <flUCtUati0n)  at  a Particular  point  depends  on 

neighboring  points  in  space  both  in  sire  and  sign.  The 

conventional  simple  linear  least  squares  method  neglects 
this  interaction,  modelling  (4.2.9)  as: 
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T = RU  + £ ' (4.2.12) 

where  £'  = S + £ is  considered  simply  as  an  error  residue 
vector. 


4 . 3 Estimation  of  the  Anomaly  Covariance  Matrix 

The  covariance  matrix  I or  the  correlation  matrix  C is  the 
quantity  which  has  to  be  evaluated  in  order  to  check  the 
validity  of  the  theoretical  predictions  of  Chapter  3.  However, 
there  are  two  difficult  problems  here  which  will  be  dis- 
cussed in  some  detail  in  the  following.  One  is  tied  to 
the  estimate  of  the  reference  value  U used  when  calculating 
the  anomalies.  If  the  "true"  value  is  U,  then 

A 

U = U + 6U  where  5U  is  an  error  which  we  do  not  know. 

We  have: 

A A 

E{ (T-RU) (T-RU) *}  = E{ (T-R(U+5U) (T-R(U+5U) ) *}  = (4.3.1) 

= E{ (S-S6U)  (S-R6U) *}  = E{ SS  * } + E{  (R6U)  (R6U) * } 

Cross-terms  E{S(R6U)*^  = E{(R6U)S*}  = 0,  because  S is 
a random  vector . The  contribution  from  the  error  covariance 
matrix  E{ (R6U) (R6U) * } in  estimating  I (eq.  4.2.10)  depend 
on  6u  as  well  as  the  coordinate  matrix  R,  i.e.,  the  fixed 
array  configuration.  Large  R may  result  in  a considerable 
error  even  if  the  elements  in  6U  are  small. 

There  are  considerable  difficulties  associated  with  the 
evaluation  of  statistical  quantities  using  a fixed  array  like 
NORSAR.  Since  the  medium  under  consideration,  the  upper  portion  of 
the  real  earth,  is  considered  random  in  space  and  not  in  time,  we 
meet  severe  sampling  problems.  Events  separated  in  time  but 
not  in  space  will  produce  the  same  response  to  the  inhomo- 
geneities, i.e.,  the  same  anomalies  will  be  observed  within 
the  limits  of  measurement  error  and  difference  due  to  fre- 
quency content  and  energy  release.  With  events  separated  in 
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space  there  is  still  a considerable  overlap  in  "sampling" 
the  inhomogeneities . By  "sampling"  overlap  or  redundancy 
sanplrng  is  meant  that  different  signals  pass  through  part  of 
the  same  structures  such  that  observations  are  not 
statistically  independent.  Statistical  independence  in  turn 
is  ruled  by  the  correlation  distance  of  the  random  medium. 

A mobile  array  would  be  able  to  overcome  most  of  these 
sampling  problems. 


In  order  to  estimate  the  covariance  matrix  ?,  10  events  with 
optimum  coverage  in  azimuth  and  velocity  were  selected.  As 
rotational  symmetry  of  the  covariance  matrix  I is  assumed, 
e e feet  of  using  the  events  with  approach  to  NORSAR  as 
epicted  in  Fig.  4.3  would  be  to  reduce  the  effect  of  strong 
local  structures  not  pertaining  to  the  randomly  inhomogeneous 
medium.  Each  event  produces  a sample  matrix  of  Z and  the 


Fig.  4.3 


Direction  of  approach  to  NORSAR  for 
the  empirical  correlation  functions 


events  used  in  estimating 
events  1-10  in  Table  1. 
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A 

average  of  10  sample  matrices  yielded  an  estimate  Z.  All 
traces  were  filtered  in  order  to  pass  through  the  energy 
at  0.7  Hz. 


In  view  of  the  estimation  problems  discussed  previously, 
it  was  decided  to  base  the  reference  U on  the  slowness  deter- 


mined from  the  NOAA  epicenter  solution.  For  logamplitudes 
the  best  reference  was  considered  to  be  the  statistical  mean 
of  all  available  sensors.  The  transformation  of  I to  6 
(Section  4.2)  needs  a few  comments.  Based  on  theoretical 
considerations,  one  would  expect  that  the  variances  of  Z to 
be  homogeneous.  This  was  not  in  fact  true,  and  we  can  think 
of  at  least  three  possible  explanations 


1)  Too  few  samples 

2)  A redundancy  sampling  effect 

3)  Local  structure  effects. 


Note  that  1)  and  2)  are  in  a certain  sense  inversely  related 
to  each  other  since  increasing  the  number  of  events  to  obtain 
better  estimates  results  in  a more  severe  redundancy  sampling 
effect.  However,  the  averaging  procedure  in  going  from  the 
132  x 132  element  matrix  C to  the  associated  correlation 
function  is  of  course  a means  to  improve  estimates  and  avoid 
redundancy  sampling,  since  different  portions  of  the  whole 
array  represent  approximately  independent  samples  of  the 
medium. 


The  correlation  functions  were  constructed  by  simply  averaging 
over  elements  in  C corresponding  to  station  separations  within 
each  [ (0. 5-N) km,  (0.5+N)km],  N=l, 2,  . . . , 100  interval.  The 
Fisher  Z-transformation  (Cramer,  1945)  for  the  correlation 
values  was  avoided  due  to  the  obvious  uncertainty  in  the  data. 
The  obtained  correlation  functions  are  depicted  in  Figs.  4.4a,b, 
together  with  the  theoretical  Chernov  (1960)  functions 
(3 . 4 . 6a,  b)  . 


SPATIAL  LA3  (torvlOI 


(b) 

Fig.  4.4  Empirical  transverse  correlation  functions  for  travel  time 
fluctuations  (a)  and  logamplitude  fluctuations  (b)  together 
with  the  theoretical  Chernov  (1960)  functions  eqs.  (3. 4. 21a, b) 
with  correlation  distance  15  km  and  wave  parameter  10. 
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4,4  Comparison  of  Observations  with  Chernov's  Theoretical 
Predictions 

The  wave  parameter  D can  be  determined  either  by  eq.  (3.4.8) 
from  an  estimate  of  oA/o4>  or  by  using  eq.  (3.4.12)  and  an 
estimate  of  the  correlation  coefficient  between  phase  and 
logamplitude.  These  two  determinations  are  independent  and 
offer,  as  point  out  by  Aki  (1973)  a severe  test  on  the 
applicability  of  Chernov's  theory  to  observed  values.  Con- 
sequently, the  quantities  a A,  a$  and  r^  were  estimated  for 
events  1-10  in  Table  1 at  0.7  Hz  and  for  events  2,  4 and  6 
at  0.9,  1.1  and  1.3  Hz.  oA  was  computed  referred  to  the  mean 
logamplitude  and  o4>  relative  to  the  best  fitting  plane  wave. 

$ is  defined  as  4>  = At*o)  where  to  is  angular  frequency  and 
At  is  the  plane  wave  residual.  Results  are  presented  in 
Table  4.  Note  that  the  number  of  sensors  used  is  seldom  132 
due  to  array  operating  discrepancies. 

Following  Aki's  procedure  to  plot  the  results  in  oA/a<J> 
versus  r^  diagram,  they  are  easily  compared  to  the  theoretical 
predicted  relationship  obtained  by  eliminating  D from  eq. 
(3.4.8)  and  eq.  (3.4.12).  Fig.  4.5  shows  how  the  NORSAR  data 
fit  into  the  picture  presented  by  Aki  for  LASA  data  using  a 
somewhat  different  measurement  procedure.  Capon's  (1974) 
average  result  based  on  273  measurements,  represented  by 
a big  black  dot,  is  also  shown.  The  lower  part  of  Fig.  4.5 
gives  the  average  of  events  2,  4 and  6 in  Table  1 for  dif- 
ferent frequencies. 

Confidence  limits  for  the  two  estimators  r^  and  ofi  o$  have 
to  be  considered  before  any  conclusions  can  be  reached. 
Introducing  the  Fisher  Z-transformation  (CramSr,  1945)  to 
the  correlation  coefficient  r$A,  it  is  shown  that  the 
variable  Z defined  as 
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Fig.  4 


.5  oA/o$  versus  r.^  diagrams  redrawn  from  Aki  (1973).  All 
text  pertains  to  his  original  figure. 

Symbol  definition 
Upper  box: 

Chernov's  theoretical  relationship 

95%  confidence  limits  for  Aki's  (1973)  observations 
(66  samples) 

.?.  Ak^-'s  (1973)  observations  for  LASA 

Confidence  limits  for  r..  estimates  in  this  stud'' 
based  on  111  samples 

Confidence  limits  for  oA/o$  estimates  in  this  study 
based  on  111  samples 

A Estimates  for  events  1-10  (0.7  Hz)  in  Table  4 
• Capon's  (1974)  average  reported  value  for  LASA  at 
0.8  Hz  (sample  size  279) 


Lower  box: 


A 

Average 

of 

events 

2, 

4 

and 

6 

at 

0.7 

Hz 

■ 

Avei age 

of 

events 

2, 

4 

and 

6 

at 

0.9 

Hz 

A 

Average 

of 

events 

2, 

4 

and 

6 

at 

1.1 

Hz 

0 

Average 

of 

events 

2, 

4 

and 

6 

at 

1.3 

Hz. 

46 


(4.4.1) 


may  be  considered  a normally  distributed  variable  with  mean 


Z 


1+r 


♦ A 


1-r 


4>A 


r 4*  A 
2 (N-l) 


and  standard  deviation 


(4.4.2) 


where  N is  the  number  of  samples  used  in  the  correlation 
coefficient  estimate. 


The  inverse  transformation  of  (4.4.1)  yields 

= tan  h (Z)  (4.4.4) 

and  the  confidence  limits  for  ^ are  obtained  by  first 
computing  the  limits  for  Z and  then  using  eq.  (4.4.4)  to 
obtain  the  limits  for  r^.  The  95%  confidence  limits  for  Z 
are  Z + 2o  and  the  limits  for  rtl  are 

r<tAl  = tan  h(Z-2oz) 

(4.4.5) 

r0A2  = tan  h(Z+2oz) 

Where  r0Al  and  r <t> A 2 denote  the  lower  and  upper  95%  confidence 
limits  for  r^A  respectively. 

The  average  cross  correlation  coefficient  r^A  estimated  for 
event  1-10  at  0.7  Hz  is  0.18  for  a sample  size  of  N=1240. 
Assume  now  that  our  HQ  hypothesis  is 


H • r 

0 ’ r$A 


0 


(4.4.6) 
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The  confidence  limits  for  are  given  by  (4.4.5)  to  be 
r4>Al  = -°*06  and  r <j>A2  = °*06 

Thus  the  probability  of  doing  an  error  when  rejecting  the 
hypothesis  (4.4.6)  is  less  than  5%  provided  the  measurement 
based  on  1240  samples  lies  outside  the  range  (-0.06,0.06). 
The  measured  value  of  r^^  was  0.18  and  is  indeed  very 
significant.  The  transverse  cross  correlation  function  for 
logamplitude  and  phase  fluctuations  was  not  developed  by 
Chernov  (1960).  However,  by  calculating  the  best  plane  wave 
phase  anomalies  and  the  logamplitude  anomalies  relative 
to  the  statistical  mean  value,  an  empirical  estimate  £ 

<J>A 

was  produced.  Assuming  as  before  rotational  symmetry,  we 
obtained  the  cross  correlation  function  depicted  in  Fig.  4.6. 


The  confidence  limits  for  the  ratio  oh/o$  are  also  interesting 
to  know.  Suppose  A and  4>  are  independent  and  normally  dis- 
tributed N (0,  o^)  and  14(0,0^)  respectively.  The  quantitites 
A and  <J>'  where 


4>'  = 


(4.4.7) 


are  still  independent  and  normally  distributed.  The  variance 
ratio 

A 2 

a i 

f 2 _ A 

t (4.4.8) 

°$' 

is  known  to  be  F-distributed  (Cramfir,  1945),  and  since 


(4.4.9) 


we  obtain 


SfiATIAL  LAG  (km  B) 


Fig.  4.6  Empirical  transverse  cross-correlation  between  phase  and 
logamplitude  fluctuations  at  NORSAR. 


f2 


A 2 

°A 


A j 

oi 


(4.4.10) 


If  the  95%  corifidence  limits  for  f2  are  denoted  y2  and  y2 
respectively,  we  have 


(4.4.11) 


which,  using  (4.4.10),  leads  to 
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(4.4.12) 


From  an  F-distribution  table  we  find  that  for  estimates  with 
NA  = N<j>  = 111  degrees  of  freedom  (cf.  Table  4),  we  have 
Y1  ~ and  y2  = 1-17.  The  confidence  intervals  for  estimates 

of  the  ratio  a is  then  within  the  curves  determined  by 


„ A 

°A  °A  ° h 

0.85  -A  < -A  < 1.17  -A 

<J>  O a<J) 

<*> 


(4.4.13) 


These  curves  are  drawn  as  thin  solid  lines  in  Fig.  4.5  and 
show  that  except  for  very  large  values  of  D (aA/o#  very  close 
to  unity)  they  lie  entirely  inside  the  confidence  limits 
for  the  correlation  coefficient  r 


As  inferred  from  Fig.  4.5,  all  NORSAR  points  are  outside  the 
confidence  region  for  Chernov's  (1960)  theoretical  predictions. 
In  order  to  see  what  effect  the  independently  determined 
NOAA  slowness  estimates  (cf.  Table  1)  have  on  the  results 
depicted  in  Fig.  4.5,  the  relevant  calculations  of 
and  r$A  were  performed  using  these  values.  Results  are  shown 
m the  upper  part  of  Fig.  4.7.  Also  the  NOAA  slowness  cor- 
rected for  a dipping  Moho  with  updip  94  deg  from  North,  dip 
6 deg  and  velocity  contrast  6 . 6/8 . 2 was  used  as  reference 
for  the  calculations.  This  particular  interface  was  shown 
by  Berteussen  (1975)  to  explain  maximum  of  the  NORSAR  time 
corrections.  Results  using  this  reference  are  shown  in  the 
lower  part  of  Fig.  4.7.  As  can  be  seen,  both  alternatives 
increase  the  estimate  and  scatter  the  data. 
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CORRELATION  COEFFICIENT  BETWEEN  AMPLITUDE  AND  PHASE 


Fig.  4.7  See  symbol  definition  of  Fig.  4.5. 

Upper  box: 

The  symbols  A give  a. /a.  versus  r. . estimates  obtained 
for  events  1-10  (Tabie  1)  at  0.7  Hz  using  slowness  estimated 
from  NOAA  epicenter  solution  and  a symmetrical  earth  model 
with  standard  velocity  distribution. 

Lower  box: 

Same  as  in  upper  box,  but  with  slowness  corrected  for 
dipping  Moho  with  updip  6 deg  in  direction  94  deg  from 
North  and  a velocity  contrast  of  6. 6/8. 2. 


Normality  in  the  distribution  of  travel  time  (phase)  fluctua- 
tions is  one  of  the  important  features  of  a wave  field  scattered 
in  an  isotropic  inhomogeneous  medium.  Fig.  4.8  shows  a histo- 
gram of  plane  wave  residuals  formed  by  the  data  from  which 
°<j>  estimates  of  Fig.  4.5  are  constructed.  It  is  quite  clear 
that  we  are  concerned  with  a skewed  distribution  which  may 
lead  to  an  over-estimation  of  values.  Berteussen  (1974) 

also  obtained  this  distribution  with  other  data. 
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Fl9’  4 6 Tab'rr.1'10"  °f  bSSL  Plari'  “ave  Nations  for  events  1-10, 


he  interesting  question  is  whether  this  skewness  has  a 
geophysical  explanation.  Local  high  velocity  bodies  within 
e crust  are  likely  to  produce  deviations  of  this  type 
n our  opinion,  such  potential  bodies  of  high  velocity  give 
oo  strong  scattering  locally  underneath  NORSAR,  thus  vie- 
wing the  conditions  of  the  Chernov  (i960)  medium  and 

LitslnVhe  V%  6StimateS  °f  F19-  4'5  >*lo«  the  confidence 
its.  Anyway,  the  hypothesis  about  random  scattering 

a ing  place  according  to  a Chernov  (1960)  type  of  medium 

has  to  be  rejected  on  the  basis  of  the  test  provided  by  the 

A/0{,  versus  r4fl  computations  in  this  study. 


e previously  referred  studies  for  LASA  data  by  Aki  (1973) 
and  Capon  (1974)  both  give  estimates  for  the  parameters  £ 
a and  L characterizing  a Chernov  (1960)  type  of  crust  and  ur 
mantle.  Aki  ,1973,  limits  his  estimation  L those  data  whic 


fell  inside  the  confidence  regions  (dashed  lines)  in  Fig.  4.5 
Capon  (1974)  used  a marginal  average  value  to  compute  his 
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parameters.  Capon  and  Berteussen  (1974)  using  NORSAR  data 
obtained  zero  correlation  between  amplitude  and  phase  fluc- 
tuations and  decided  that  the  scattering  beneath  NORSAR  was 
too  strong  for  the  Chernov  type  of  medium  to  be  valid  in  the 
frequency  range  0.6  to  2.0  Hz.  As  mentioned  previously, 
these  works  were  based  on  spectral  data  and  somewhat  dif- 
ferent measurement  procedures.  Recently  Berteussen  et  al 
(1975)  estimated  Chernov  (1960)  parameters  from  the  data 
prepared  for  this  study,  using  the  method  described  by 
Chris toff ers son  (1975).  However,  this  method  does  not  include 
the  o^/Oq  versus  r$A  test  on  the  applicability  of  Chernov's 
theory,  and  is  rather  based  on  the  hypothesis  that  Chernov's 
theory  is  valid.  The  hypothesis  could  not  be  rejected  from 
the  parameter  values  obtained. 

The  phase  fluctuations  beneath  NORSAR  are  obviously  too 
strong  to  accept  that  the  crust  and  upper  mantle  beneath 
NORSAR  act  like  a Chernov  medium.  On  the  other  hand,  the 
"abnormal"  phase  fluctuations  do  not  seem  to  violate  the 
predicted  positive  cross  correlation  r$A>  As  inferred  from 
Table  2 r^  = 0.18  in  average  for  event  1-10  in  Table  1, 
which  is  a reasonable  value.  Assuming  an  uncertainty  of 
about  one  kilometer,  we  read  from  Fig.  4.4  that  the  cor- 
relation distance  is  about  15  km.  With  a mean  velocity 
aQ=8  km/sec  in  the  crust  and  upper  mantle,  corresponding  to 
k=0.55  rad/km  at  0.7  Hz,  we  are  able  to  determine  the  depth 
L from  the  definition  of  D (eq.  3.4.7).  As  we  can  see,  the 
Fresnel-approximation  of  ka  being  much  larger  than  1.0  is 
supported  by  a value: 

ka  c*  8 2 

From  eq.  (3.4.12),  we  obtain: 

D = 15  >>  1 


(rad) 

0.48  1.06 

0.56  1.15 

0.39  0.93 

0.33  1.22 

0.42  0.70 

0.34  0.95 

0.33  1.08 

0.40  1.00 


ua/uv  r . 
(rad"  ) <r>A 


0.23 
0.  19 
0.23 
0.  31 
0.  15 
0. 12 
0.  31 
0.  22 


Number  of 
samples 


TABLE  2 

Root  mean  square  fluctuations  in  phase  <t  and  logamplitude  A 

/e  9iIIh  t0gfher  "lth  the  correlation  between  these  fluctuations 

f ° Sdmples  corresponds  to  number  of  sensor  recordings 
available  for  a particular  event.  y 


mm 
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The  depth  of  the  inhomogeneous  medium  is  then  fixed  to 
L c*  464  km  > a ai  15  km 

The  mean  square  fluctuations  n2  of  the  refractive  index  can 
be  determined  either  from  an  estimate  of  A2  and  eq.  ( 3 . 4 . 6 . b) 
or  P and  eq.  (3. 4. 6. a) . 

The  former  quantity  gives  (see  Table  2) 

n2  = 0.000095  <<  1 (4.4.1 

The  Born-approximation  is  valid  for  small  fractional  energy 
losses  AI/I  where: 

AI/I  = /?  n2  k2a  L(l-e"k  a ) (4.4.1 

Computing  this  quantity  using  logamplitude  fluctuation 
yields 


AI/I  = 0.35 


(4.4.14) 


which  is  on  the  edge  of  an  acceptable  value.  On  the  other 
hand,  using  phase  fluctuations  for  determining  the  refractive 
index  fluctuations  we  get: 


n2  = 0.00046 


(4.4.15) 


and 


* i'74  (4.4.16) 

Conclusively,  we  infer  that  the  theoretical  formulae  lead  to 
a contradiction  when  phase  information  are  used  in  the  esti- 
mation of  refractive  index  fluctuations  and  fractional  energy 
loss  (cf.  eq.  (4.4.16)).  If  we  consider  the  less  restrictive 
Rytov-approximation  done  to  obtain  the  formulae  in  Chapter  3, 
requiring  that  the  energy  lost  per  wave  length  be  small,  we 
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obtain  values  satisfying  our  assumptions.  Finally,  we  observe 
from  the  lower  part  of  Fig.  4.5  that  increasing  frequency 
is  likely  to  correspond  to  lower  values  of  D in  agreement 
with  the  formula: 

2 LX 

D = 4L/ka2  = - — {X  denotes  wave  length)  (4.4.17) 

a2 

On  the  other  hand,  the  applicability  of  the  Chernov  theory 
is  not  improved  for  higher  frequencies. 
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5-  THE  stochastic  term  in  the  refined  model  for  time  and 
LOG AMPLITUDE 

5 • 1 Generalized  Least  Squares  Estimation  and  Prediction 

The  transver:  ? autocorrelation  functions  shown  in  Figs. 
4.4a,b  give  striking  evidence  for  the  stochastic  term  intro- 
duced in  paragraph  4.2.  Conventional  modelling  according  to 
the  simple  least  squares  model  represented  by  eq.  (4.2.12) 
presupposes : 


E { € ' } = 0 


(5.1.1) 


E { 6 ' 6 ' * } = o 2 , I 
e 


(5.1.2) 


(I  is  the  identity  matrix),  and  the  essential  constraint  is 


n 

minimum  (7  e? ) 

i = l 1 


(5.1.3) 


The  problem  with  correlated  residuals  can  be  circumvented 
by  linear-transforming  observations  using  the  covariance 
matrix  I = E{SS*}.  Disregarding  abnormal  cases,  I is  a posi- 
tive definite  matrix.  Then  it  is  invertible  and  expressible 
as  (Ferguson,  1967) 


* ' AA  (5.1.4) 

where  X is  also  a positive  definite  invertible  n x n matrix. 

Including  the  negligible  error  term  €'  in  the  S-term  and  trans- 
forming our  variables  through  X-1  yields  for  eq . (4.2  . 9 ): 

X T = X RU  + X ^ S it  i 


Now,  since 
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E{ (X-1T-E{X“1T})  (X-1T-E{X_1T} ) *}  = E{ (X-1S) (X_1S) * } 

X-1E{SS*}X-1  = X-1EX_1  = X_1XXX“1  = I (5.1.6) 

our  observations  have  been  transformed  into  uncorrelated  data 
with  qualifications  satisfying  eq.  (5.1.1)  and  5.1.2).  The 
simple  least  squares  technique  can  now  be  applied  to  the 
variables  requiring  that  we  perform  minimization  of: 

(X_1S) * (X_1S)  = S*I_1S  (5.1.7) 

The  covariance  or  correlation  matrices  of  our  data  provide 
information  which  can  be  taken  advantage  of  in  order  to 
predict  values  at  certain  points  in  space.  At  this  stage  it 
is  convenient  to  introduce  the  combined  procedure  of  gen- 
eralized least  squares  estimation  and  prediction  as  actually 
done  in  Dahle  et  al  (1975).  Assume  coefficients  for  the 
linear  trend  have  to  be  estimated  from  n basic  observations, 
and  that  the  observed  field  is  to  be  predicted  at  m other 
points.  Including  as  before  the  measurement  error  term  £'  in 
the  stochastic  term  S,  we  obtain,  using  matrix  designations 
and  the  concept  of  partitioned  matrices: 


T = RU  + BS 


where 


R 


1 


lx 


r 

nx 


ny 


(5.1.8) 


(5.1.9) 


(5.1.10) 


1 


r 


- 5«  - 


U = U 


(5.1.11) 


B ” ^nn'  0nm  1 = n TOWS 

nn  nm 


1 0 ...  0 0 


0 1 


1 0 

0 10 


(5.1.12) 


n columns  m columns 


(5.1.13) 


And  provided  that  the  covariances  matrix 


Z = E { SS  * } = 


(5.1.14) 


is  known,  the  estimates  0 and  £ are  given  by: 


0 = m‘1r*n'1t 


(5.1.15) 


§ = ZB*N_1(T-rC) 


(5.1.16) 


where 


N = BZB* 


(5.1.17) 


= R*N_1R 


(5.1.18) 


f 
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Then  we  obtain 


N = [ I ,0  ] 

nn  nm 


nn  nm 


mn  mm 


nn 


nm 


= I (5.1.19) 

nn 


IB' 


d I 1 

nn  nm 

'l  ‘1 
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nn 

nn 
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nm 

mn 

mn  mm 

L J 

-* 

(5.1.20) 


We  have: 
0 = 


(R*I  Ir)*1  R*I~^T 
nn  nn 


(5.1.21) 


which  is  the  generalized  least  squares  formula  for  estimating 
the  linear  trend  in  the  n basic  observation  points  (Johnston, 
1963) . We  also  see  that 


$ = 


1 

— 

I 

n 

nn 

s 

I 

m 

mn 

I i (T-rO) 
nn 


This  gives 


3 = t-rO 

n 


as  we  would  have  expected,  and 

§ = I I-1  (T-rO) 

m mn  nn 


(5.1.22) 


(5.1.23) 


(5.1.24) 


which  is  the  prediction  formula  for  the  m points  to  be 
predicted . 


Fig.  5.1  Basic  points  (crossed)  and  points  to  be  predicted 
(open) . 


The  general  all-to-all  dependence  which  exists  in  the  travel 
time  and  amplitude  fluctuations  of  a network  of  stations 
like  NORSAR  deserves  considerable  attention.  The  estimation 
and  prediction  theory  for  dependent  variables  provides  a con- 
venient tool  to  study  the  significance  of  correlated  data, 
being  stochastic  fluctuations  or  not.  The  special  procedure 
undertaken  was  to  randomly  select  half  of  the  54  stations 
comprising  subarrays  01A,  01B-07B  and  06C  as  basic  observa- 
tion points,  trying  to  predict  the  observed  values  at  the 
remaining  ones  by  estimating  eq.  (5.1.21)  and  eq.  (5.1.24). 
Basic  points  are  marked  with  crosses  in  Fig.  5.1,  while  the 
open  symbols  are  points  where  travel  time  (or  logamplitude) 
has  to  be  predicted.  The  quantity  used  as  a goodness  of  fit 
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of  the  data  to  the  stochastic  model  represented  by  eq . (4.2.9) 
is  sum  of  squared  differences  between  observed  and  predicted 
values,  a quantity  we  shall  prefer  to  denote  prediction  error 
variance . 

In  order  to  obtain  a meaningful  comparison,  the  same  quantity 
using  the  simple  least  squared  determined  plane  wave  solution 
is  computed.  This  is  nothing  else  than  the  conventional  sum 
of  squared  plane  wave  residuals  when  we  speak  of  travel  time. 
The  ratio  of  the  prediction  error  variances  in  per  cent  is 
a measure  of  the  applicability  of  the  stochastic  model,  in  such 
a way  that  a hundred  per  cent  means  no  improvement  relative 
to  simple  least  squares,  while  zero  per  cent  means  complete 
agreement  with  the  stochastic  model. 

As  inferred  from  Chapter  3,  practical  application  of 
eq.  (5.1.21)  and  (5.1.23)  involves  implementation  of  the 
matrix  I or  its  related  matrix  C.  The  latter  is  easily 
obtained  from  eqs.  (3.4.21a,b)  or  alternatively  from  the 
empirical  correlation  functions. 

Since  the  Chernov  functions  depend  essentially  on  correlation 
distance  and  the  wave  parameter  D,  an  iterative  computer  pro- 
gram v as  constructed  in  order  to  determine  the  pair  (a,D) 
which  gives  the  minimum  prediction  error  variance.  Due  to 
the  very  time-consuming  matrix  inversions  involved  it  was 
considered  unrealistic  to  try  to  obtain  a minimum  prediction 
error  variance  based  on  some  convergence  criterion.  Instead, 
after  some  tentative  runs  on  the  computer,  it  was  decided 
to  construct  a grid  of  points  (a,D)  for  reasonable  values 
of  these  quantities  and  compute  the  relative  prediction  error 
variance  for  these  points.  The  distance  between  grid  points 
was  chosen  as  3 km  for  correlation  distance  and  10  for  wave 
parameter.  Based  on  the  calculated  grid,  a contour  plot  was 
drawn  on  the  computer  by  linear  interpolation  contouring 
the  lines  of  constant  relative  prediction  error  variance. 
Results  for  three  events  are  shown  in  Fig.  5 . 2a , b , c , d , e , f 
for  different  frequencies. 
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The  minimum  prediction  error  variance  reflects  the  form  of 
the  symmetric  correlation  function  best  fitting  the  data.  As 
seen  from  the  contour  plots  of  Fig.  5.2,  there  is  no  mimrmum 
P n in  the  range  of  (a,D)  values  examined.  We  observe 
however,  that  there  is  a minimum  trough  running  parallel 
to  the  D-axis,  leaving  this  value  undetermined  but  giving 
quite  well-defined  correlation  distance.  Noteworthy,  up 

to  90%  of  the  plane  wave  residual  variances  can  be  explained 
as  a stochastic  "scattering"  effect. 

Pig.  5.3  shows  the  distribution  of  possible  combinations 
versus  horizontal  station  separation  for  the  54  instruments 
used  in  the  prediction  procedure.  Since  there  is  a definite 
minimum  of  information  around  10  km,  we  could  attribute 
the  low  correlation  distance  to  a bias  effect  caused  by 
nonuniformity  in  this  distribution. 


SPACING  IHfl) 


Fig.  5.3 


Numbpr  of  possible  instrument  combinations 
for  subarrays  01A,  0lBf...,07B,  OfoC. 


versus  spacing 
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Among  those  correlation  functions  tried,  which  were  of  the 
type  exponentials  multiplied  by  polynomials,  the  Chernov 
type  gave  the  best  over-all  performance.  The  empirical 
functions  depicted  in  Fig.  4.4  were  also  smoothed  and  intro- 
duced in  the  prediction  procedure,  but  with  limited  success. 
On  the  other  hand,  these  functions  should  not  be  used  unless 
all  instruments  were  included  in  the  prediction  experiment 
since  they  are  estimated  for  the  whole  array  area. 


5,2  Seismic  Velocity  - A Sample  from  Some  Statistical 
Distribution 

Measurement  of  a logarithmic  amplitude  should  be  associated 
with  drawing  a sample  from  a Gaussian  distribution.  This  fact 
was  mentioned  already  in  the  introduction,  and  is  a valuable 
assertion  made  use  of  in  analyzing  seismic  data  (Ringdal, 
1974).  In  view  of  seismic  wave  scattering,  depending  on 
what  part  of  the  wave  front  we  consider  across  a network 
of  stations  like  NORSAR,  the  wave-front  normal  points  in 
different  directions.  Thus  measuring  velocity  and  azimuth 
(slowness)  using  stations  with  aperture  of  the  order  na2 
(a  is  correlation  distance) , corresponds  to  estimation  of  a 
’local'  wave  front  normal  and  its  associated  velocity. 

Mack  (1972)  introduces  a wavenumber  distribution  to  explain 
seismic  signal  coherency  fall-off  with  sensor  separation. 

His  approach  agrees  well  with  the  idea  of  non-planar  wave 
fronts  caused  by  small  irregularities  randomly  distributed 
in  the  structure  considered.  The  conventional  velocity  esti- 
mation technique  applies  to  eq . (4.2.12)  minimizing  the 
residuals  €'  assumed  to  be  uncorrelated  in  space.  In  view 
of  the  previous  sections,  eq . (4.2.9)  represents  a far  more 
realistic  model  for  the  teleseismic  P travel  times.  The 
difference  between  the  two  methods  is  discussed  in  Dahle 
et  al  (1975)  where  it  is  pointed  out  that  the  choice  between 
them  is  ruled  by  sensor  separation  and  array  aperture  related 
to  the  correlation  distance.  In  statistical  terms,  generalized 
least  squares  is  a more  'effective'  model  when  sensor  separa- 
tion is  small  and  the  array  aperture  not  too  large  compared 


to  the  correlation  distance.  Grouping  the  subarrays  of 
NORSAR  into  seven  subsets  according  to  the  scheme 

Subset  Subarrays 

* OlA, OlB, 07B 

OlC , 02C, 14C 
f11  02B,03C,04C 

,,V  03B , 05C , 06C 

„T  04B , 07C , 08C 

05B, 09C , IOC 
VI1  06B, 11C, 12C 

provides  the  opportunity  to  obtain  seven  measurements  of 
s owness  for  the  same  event.  Slownesses  were  measured  both 

“S1"9 M-2-12)'  whi<*  "=  «11  Model  I,  and  using  eq. 
4.2.9)  denoted  Model  II.  I„  the  latter  case  a flxed  CQrrel 

ion  matrix  was  used,  calculated  from  eq.  (3.4.21a)  with 

m and  0=10.  This  correlation  matrix  was  considered  to 
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Fig.  5.5  Same  as  Fig.  5.4  for  Model  n. 


be  an  optimum  choice  based  on  the  prediction  experiment 
results  for  the  0.7  Hz  data  of  Fig.  5.2.  The  results  ob- 
tained are  shown  in  Fig.  5.4  for  Model  I and  Fig.  3.5  for 
Model  II.  The  effectiveness  of  the  two  methods  was  compare 
using  the  area  of  the  95%  confidence  ellipse  calculated 
from  the  estimated  slownesses  (Dahle  et  al,  1975).  The  im- 
provement of  Model  II  relative  to  Model  I is  not  easily 
seen  by  comparing  Fig.  5.4  and  5.5  visually,  but  the  numer 
comparisons  were  conclusive  . The  U-space  solutions 
presented  in  Figs.  5.4  and  5.5  reflect  the  complexity 
in  the  P-wave  travel  times  across  a seismic  network  like 
NORSAR,  and  show  that  application  of  small  aperture  seismi. 
array  measurements  provides  unreliable  slowness  estimates 


so  that  a spatial  sampling  by  a number  of  such  arrays  seems 
necessary  in  order  to  estimate  an  averaqe  wave  front.  The 
average  slowness  (stars)  of  Figs.  5.4  and  5.5  compared  to 
the  NOAA  symmetrical  earth  solution  (square)  show  the  same 
picture  as  NORSAR's  region  correction  vectors  in  Fig.  5.6, 
while  individual  estimates  (numbers)  provide  the  possibility 
of  different  interpretations. 


Fig.  5.6  Location  calibration  vectors  plotted  in  slowness  space  The 
tail  of  the  arrow  represents  the  observed  point,  while  the 
head  represents  the  NOAA  solution. 
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6 . CONCLUSIONS 

The  main  purpose  of  this  study  has  been  to  demonstrate  the 
salient  features  of  the  fluctuations  in  seismological  para- 
meters as  experienced  at  NORSAR,  and  to  introduce  an  alterna- 
tive to  the  homogeneous  earth  model.  There  is  limited  success 
with  regard  to  the  proposed  scattering  hypothesis,"  however, 
in  our  view  the  scattering  approach  deserves  considerable 
attention.  Such  a conclusion  is  justified  by  the  fact  that 
the  seismic  wave  field  exhibits,  to  a satisfactory  degree, 
some  of  the  statistical  qualities  predicted  by  the  theory. 

The  assertion  of  forward  scattering  and  mode  conversion 
being  basically  P-P  turns  out  to  be  an  important  point. 

As  indicated  in  Section  4,  there  is  reason  to  suspect  the 
crust  and  upper  mantle  locally  being  inhomogeneous  in  a way 
not  pertaining  to  Chernov's  (1960)  theory  for  weak  inhomo- 
geneities equally  distributed  throughout  the  medium.  Stronger 
inhomogeneities  giving  rise  to  back-scattering,  conversion 
to  modes  other  than  P and  higher  order  scattering  (scattered- 
scattered  waves)  are  possible  explanations  of  the  discrepan- 
cies between  theory  and  experience. 

Substantial  support  to  the  scattering  hypothesis  is  provided 
by  the  correlation  functions  depicted  in  Figs.  4.4  and  4.6. 
These  are  the  features  leading  to  the  space  predictive 
procedure  outlined  in  Chapter  5,  which  also  has  been  used  in 
interpolation  and  extrapolation  of  point  gravity  measure- 
ments (Heiskanen  and  Moritz,  1967).  The  prediction  procedure 
is  important  for  two  reasons.  Firstly,  because  it  shows  the 
anomaly  percentage  explainable  by  correlated  stochastic 
terms,  and  secondly,  it  provides  a means  of  finding  the 
optimum  correlation  matrix  for  the  data  (c.f.  Fig.  5 ?). 

The  realization  of  randomly  corrugated  seismic  wave  fronts 
is  undoubtedly  important  in  small  aperture  array  seismology. 

We  have,  in  our  opinion  successfully,  related  the  fluctuations 
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in  body  wave  travel  time,  local  wave  front  normal  and 
amplitude  to  wave  scattering  in  an  inhomogeneous  crust  and 
upper  mantle.  A quantitative  description  of  the  inhomogeneous 
medium  is  inconclusive  for  the  moment,  but  our  results  show 
that  attention  must  be  paid  to  scattering  effects  when  ob- 
served times  or  amplitudes  are  inverted  into  structural 
models . 

Finally,  as  a comment  to  the  large  amplitude  fluctuation 
across  a seismic  network  like  NORSAR,  we  see  from  Table  2 
that  the  average  standard  deviation  of  logampli tudes  is 
0.4.  This  number  gives  rise  to  a standard  deviation  in 
magnitude  measurements  (assuming  constant  period)  as  for- 
warded by  Veith  and  Clawson  (1972)  for  a world  wide 
network  of  stations  reporting  to  USCGS. 
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